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Abstract 

Recent experiments and numerical simulations have shown that certain types of microorganisms 
“reflect” off of a flat surface at a critical angle of departure, independent of the angle of incidence. 
The nature of the reflection may be active (cell and flagellar contact with the surface) or passive 
(hydrodynamic) interactions. We explore the billiard-like motion of a body with this empirical 
reflection law inside a regular polygon and show that the dynamics can settle on a stable periodic 
orbit or can be chaotic, depending on the swimmer’s departure angle and the domain geometry. 
The dynamics are often found to be robust to the introduction of weak random fluctuations. The 
Lyapunov exponent of swimmer trajectories can be positive or negative, can have extremal values, 
and can have discontinuities depending on the degree of the polygon. A passive sorting device is 
proposed that traps swimmers of different departure angles into separate bins. We also study the 
external problem of a microorganism swimming in a patterned environment of square obstacles, 
where the departure angle dictates the possibility of trapping or diffusive trajectories. 
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1. Introduction 

Microorganisms often follow unexpected trajectories when swimming near surfaces. Organisms 
such as spermatozoa and E. coli are attracted to surfaces and may accumulate there due to a 
combination of hydrodynamic and steric effects [liaEliaElElE!- This has biological and even 
medical implications, as attraction and accumulation may lead to the development of biofilms 
[EiEiiini, and infection of medically-implanted surfaces m- The orientations of swimming bodies 
in such an environment depend on the geometry of the swimmers, their mechanism of propulsion, 
and their interactions [I2iiisiii2iisiiiniii7i[ii[i3[7i[2n]. In addition to being attracted to walls, 
E. coli cells tend to swim in large circular trajectories due to hydrodynamic interactions with the 
surface [2T]; notably, the orientation of the circular trajectories are reversed for swimming near a 
fluid-air interface [221123] , which can be rationalized by considering the two different image flows 
near no-slip and shear-free surfaces [MHZ]. 

The unusual behaviors of microorganisms near surfaces has motivated the development of many 
intriguing engineering applications. Wedge-shaped obstacles have been designed to passively sort 
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E. coli cells into areas of differing concentrations |24j . and E. coli cells have been used to drive 
micro-devices powered by asymmetric gears [251 [261 Ea EH]. Sorting and rectification devices that 
exploit the interactions of microorganisms and asymmetric surfaces (including funnels and gears) 
have also been designed and studied [231 ESEHIESIISI] . In some cases, steric collisions or near-held 
lubrication forces may dominate long-range hydrodynamic effects 

Swimming trajectories are naturally more intricate in these complex environments. For in¬ 
stance, Takagi et al. showed that microswimmers in a held of passive colloidal beads display a 
billiard-like motion between colloids, intermittent periods of entrapped, orbiting states near sin¬ 
gle colloids, and randomized escape behavior [33|. Brown et al. extended this work to swimming 
through a “colloidal crystal,” where a synthetic swimmer hops from colloid to colloid with a trap¬ 
ping time that depends on fuel concentration, while E. coli trajectories are rectihed into long, 
straight runs [35| . Coupling topography to propulsive mechanism has also been exploited to guide 
active Janus particles |36j . A far-held hydrodynamic model predicts attraction, entrapment, and 
scattering of particles near spherical colloids 133, though the nature of entrapment is generally 
specihc to the microorganism or active particle [HTlEHlIMlIlQlSIlSZliSI . E. coli, for instance, 
can be trapped by convex walls in part because the cell body pitches down toward the surface in 
equilibrium mini; moreover, the tumbling behavior of E. coli is suppressed near surfaces due to 
increased hydrodynamic resistance |44[ I45| . For a review of the dynamics of active particles in 
complex environments see |46| . 

Another recent experiment showed very clearly that the wall interaction may depend sensitively 
upon the microorganism geometry: Kantsler et al. found that Chlamydomonas algae cells scatter 
from a hat wall due to contact between its hagella and the surface (Fig. [^), so that the interaction 
is highly dependent on the cell body shape and hagellar lengths |37|. With each stroke, one or both 
hagella are in contact with the surface and the cell body rotates as a consequence of the resultant 
torque. Rotation continues until there is no further hagellar contact with the surface, at which 
point the cell swims back into the bulk huid. As reproduced in Fig. the study included an 
examination of mutant cells with hagella of different length, and cells were found to depart from the 
surface on average at a critical rehection angle predicted by a simple geometrical consideration; the 
departure angle of wild-type Chlamydomonas was reported to have a narrow peak around 9c = 16°. 
The nature of the interaction suggests that the departure angle relative to the wall is independent 
of the angle of incidence, distinguishing this interaction from that of a classical billiard-like specular 
rehection. This angular rectihcation was also observed by Spagnolie &: Lauga in an analytical study 
of a model “potential how squirmer” (a model microorganism with a prescribed tangential surface 
velocity) swimming near a surface [7|. 

Motivated by these behaviors of microorganisms and active particles near walls, we explore the 
billiard-like motion of a swimming body which, upon wall impact, swims away with a hxed angle 
of departure independent of the angle of incidence. While classical billiard dynamics — in which a 
particle undergoes symmetric (specular) rehections with each wall collision — is a well-developed 
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Figure 1: (a) A Chlamydomonas reinhardtii cell approaches a flat wall. Flagellar interactions with the surface lead 
to cell rotation until the cell swims back into the bulk at a geometry-specific departure angle, (b) The departure 
angle depends on the flagellar length, supporting the geometric view of the wall interaction. Reproduced from m 
with permission. 


problem in mathematics [H], the microorganism billiard motion of our investigation has not to our 
knowledge been studied even in simple domains. A similar type of non-specular reflection law also 
appears in spiral waves in bounded media |1U], cavity solitons m. optical microcavities m. and 
even in internal gravity waves in stratified fluids [52l |53l [5lj . Techniques used to study the classical 
billiard system, such as periodic reflections of the fundamental domain to tile a surface |55j . hinge 
upon symmetry not present in our system. In the mathematics literature, nonstandard reflection 
laws have recently been studied [SniEZlEH], but the focus has been mostly on contracting laws. 
The slap map [58] applies to our setting for the very limited case of perpendicular reflection. 

We begin by studying two-dimensional microorganism billiard dynamics inside a regular poly¬ 
gon of arbitrary degree, where a one-dimensional return map is used to characterize the possible 
dynamics into stable periodic trajectories or unstable, chaotic orbits. We show that the dynamics 
are robust to the introduction of weak random fluctuations. The Lyapunov exponent describing 
the dynamics of an ensemble of swimmer trajectories can be positive or negative, can have ex¬ 
tremal values, and can have discontinuities depending on the degree of the polygon. Using the 
results for the stable periodic trajectories we propose a passive sorting device to trap swimmers of 
two different departure angles into two separate bins. Finally, a similar method of investigation is 
applied to the external problem, wherein the body swims through a periodic array of square ob¬ 
stacles. The results may contribute to our understanding of swimming microorganisms in porous 
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environments, and suggest potential applications in entrapment and passive hydrodynamic sorting, 
and bioremediation. 

The paper is organized as follows. In ^ we study microorganism billiard dynamics inside a 
regular polygon, where consecutive and non-consecutive wall impacts are considered. An explicit 
formula for the stable periodic trajectory is derived, and is shown to be robust to the addition of 
weak random fluctuations. The invariant measure and the Lyapunov exponent for the trajectories 
are also discussed. A sorting device is then proposed in ^ Finally, in ^we turn our attention to 
the external problem, billiard-like swimming in a periodic field of square obstacles, where lessons 
from the interior problem are instructive but the particle trajectories are more widely varied. We 
conclude with a discussion in ^ 
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Figure 2: (a) Schematic representation of microorganism billiard wall interactions. Regardless of the incident angle, 
8i, the swimmer slides along the wall a distance S and then departs with angle 6c. (b) Microorganism billiard inside 
a regular hexagon with interior angle 9p = 27r/3. The nth point of impact, Xn, is measured from the vertex just 
behind the swimmer. 
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2. Microorganism billiards inside a regular polygon 

2.1. Adjacent wall reflections lead to stable periodic orbits 

We begin by studying the deterministic microorganism billiard inside a regular polygon with 
N sides, in the case where the body can only move from each wall to an adjacent wall. Figure 
illustrates the collision interaction with any wall; regardless of the angle of incidence, 9i, the body 
departs from the surface at a critical departure angle, 9c G (0, vr/2). We include the possibility that 
the body translates a distance 6 between impact and departure. The constraint that the swimmer 
only impacts an adjacent wall can be written in terms of the departure angle; namely, with the 
interior angle of the polygon given by 9p = {N — 2 )tt/N, we require 9c G (0, n/N). 

To understand the nature of the dynamics we take advantage of the symmetry in the system 
and study a map of [0,1] to itself, writing Xn for the nth point of impact, measured from the vertex 
just behind the swimming cell (see Fig. Ep)- (This is called the reduced billiard map [58].) Given 
Xn, the subsequent point of impact is found geometrically. Assuming for the moment that <5 = 0, 
we find 


Xn+l = f{Xn) = l3{l - Xn), (la) 

fl = sin(6'c)/ sin(27r/Ar - 9c), (lb) 


where /3 < 1 in the case of adjacent wall impacts. Writing Xn = /3(1 — Xn-i), and so on, results in 
an explicit expression for the nth position. 


Xn = i-fl)''xo - 

i=l 

1 _ (- 8 )^ 

= i-flrxo+p—^^8^. ( 2 ) 

Since fl < 1, the dynamics lose memory of the initial position xq exponentially fast in the number 
of impacts and the body settles to a periodic orbit with 


x* = lim Xn 

n^oo 


i + fl’ 


(3) 


which is always less than 1/2. Consider a uniform collection of swimmers departing from the entire 
length of one wall. By Eq. ( |la| ), the swimmers travel as a “beam” through the domain until colliding 
with the adjacent wall in a smaller region of length /3, so that fl may be viewed as a stretch (or 
compression) factor. In the present setting the beam of swimmers is focused into a smaller region 
at every iteration, so that the trajectory rapidly converges to an asymptotically stable periodic 
orbit. Examples of this situation are shown in Fig. [^<feb, with 9c = 30° and 9c = 35.3°, both less 
than 36° and so resulting solely in adjacent wall reflections. The stable periodic orbits are, as they 
must be, inscribed pentagons, and the rate of convergence depends on the stretch factor (3. In the 
second case fl is very nearly one, and the convergence is slow. 


5 





Figure 3: (a) The trajectory with departure angle 6c = 30° and initial point xq = 0.1 inside a pentagon rapidly 
converges to the stable periodic orbit. The trajectory is shown in red after the first 15 reflections, (b) The trajectory 
with departure angle 6c = 35.3° and initial point xo = 0.1 converges more slowly to the stable periodic orbit. The 
trajectory is shown in red after the first 150 reflections. 


However, if we allow for the boundary case 9c = then the full length of the initial wall 

maps to the full length of the adjacent wall. Correspondingly we have /? = 1, and instead of 
the trajectory settling to one stable periodic orbit we now find that every initial point xq resides 
on a neutrally stable periodic orbit. For motion inside a square, for instance, the special angle is 
9c = 7 r/ 4 , and any inscribed rectangle represents a neutrally stable trajectory. An example is shown 
for motion inside a pentagon with xq = 0.1 in Fig. & as a dashed line. Generally, an even number 
of walls leads to periodicity with a fundamental period (i.e., with N — 1 reflections), while an odd 
number of walls leads to periodicity with twice the fundamental period (with 2N — 1 reflections). 
The initial position xq = 1/2 is a special case, for which the swimmer returns to the initial point 
after — 1 reflections for any N. 

Returning to the case that <5 7 ^ 0, so that the body slides along the wall a hxed distance 5 
after each wall collision before departing, the mapping from the departure point Xn to the next 
departure point Xn+i is easily modified to Xn+i = /?(! — Xn) + S (assuming that 6 is small enough 
that the swimmer does not slide past a vertex, Xn+i < 1 ), resulting in the trajectory 

= + + (4) 

Hence, the motion still settles to an asymptotically stable periodic orbit in this case, with a shift 
in the fixed point of departure to x* = (5 + /?)/(! + /?). More generally, the sliding distance 
6 might depend on the angle of incidence, as for model potential flow squirmers [7] and model 
“pusher” particles m- However, under the assumption of adjacent wall reflections, since the angle 
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Figure 4: When 6c = kn/N, each initial point xo is part of a neutrally stable periodic orbit. Shown are the (a) 
pentagonal, (b) hexagonal, and (c) octagonal boundary geometries, and trajectories for k — 1 (dashed lines), k = 2 
(solid lines), and in the last case, k = 3 (dotted line), all with xq = 0.25. Arrows indicate the swimming direction. 


of incidence is independent of the departure position Xn, then 5 is in fact constant, and the previous 
theory applies. 

2.2. Non-adjacent wall reflections: periodic and chaotic dynamics 

We will return to the role of random fluctuations in the dynamics, but first we will show that 
strong disorder in the system can be generated by the dynamical system alone once the constraint 
of purely adjacent wall reflections is relaxed. We now allow that the swimmer’s next wall of impact 
may depend not only on 8c but also on its position Xn- Fortunately, from the perspective of 
mathematical tractability, for a given departure angle a swimmer departing from any one wall can 
only possibly reach at most two other walls in the deterministic setting (which can be understood 
geometrically by inscribing isosceles trapezoids inside the regular polygon). For this study we set 
the sliding distance 5 to zero, and we will only consider departure angles with 9c < t:{2. 

For one of a set of special departure angles, Oc = Trk/N, where k G {1, 2,..., N — 2}, the beam 
of swimmers collides with only one other wall. For k < {N — l)/2, the map is given trivially by 
Xn+i = 1 — Xn. Just as was seen for the adjacent wall case {k = 1), there are an infinite number 
of neutrally stable periodic trajectories, since Xn +2 = Xn. Let pk = mN, where p and m are the 
smallest integers such that this relation holds. Then any initial point xq corresponds to a neutrally 
stable orbit with a period (measured in the number of hits) equal to p if p is even or xq = 1/2, and 
period 2p if p is odd and xq 7^ 1/2. 

Sample trajectories beginning at xq = 0.25 inside the N = 5, N = 6, and N = 8 polygonal 
domains are shown in Fig. |^-c, for k = 1 (dashed line trajectories), k = 2 (solid line trajectories), 
and in the last case, k = 3 (the dotted line trajectory). For N = 5, p = 5 for both k = 1 and 
k = 2, so that a full period requires 10 reflections. For N = 6, p = 6 for both k = 1 and k = 2, 
so that a full period requires 6 reflections. Finally, for N = 8, p = 8 for k = 1 and k = 3, while 
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p = 4 for k = 2. In this last case the neutrally stable periodic rectangle inscribed inside a square 
appears, as the dynamics do not explore the regions where the octagonal boundary differs from a 
square boundary, and we are back to the adjacent wall case. 



Figure 5: (Left) When N is even, the swimming direction at the next reflection requires an orientation reversal 
for the map when 9c G {{N — 2)n/2N,n/2), as indicated for the case N = 8. Two sample trajectories leaving 
from the bottom wall with 9c = 76° are shown in red. (Right) When N is odd, a reversal is required when 
9c G {{N — 3)'k/2N, {N — l)7r/2A'^) and when 9c G {{N — l)7r/2A'’, 7r/2), as indicated for the case N = 7. Two sample 
trajectories leaving from the bottom wall with 9c = 56° are shown in red. 


Meanwhile, an additional subtlety must be taken into account for a small range of departure 
angles, specifically those for which the swimmer, as it approaches the next wall of impact, has 
changed its orientation relative to the mapping scheme: the direction of reflection must always 
be towards increasing x. When N is even, the swimming direction at the next reflection requires 
an orientation reversal for the map when 9c € {{N — 2)7r/2N,7r/2). When N is odd, a reversal is 
required when 6c G {{N — 3)7r/2A^, [N — 1)'k/2N) and when 9c G {{N — l)7r/2Af, 7r/2). Figure]^ 
shows the cases N = 7 and N = 8. However, every wall may still be identified with every other 
wall if the direction of increasing x is always the direction of swimming, also indicated in Fig. 

For departure angles in between these special values, a beam of swimmers extending from the 
full length of the original wall is bisected and collides with two adjacent walls. The map from 
Xn to Xn +1 is piecewise linear with at most two distinct domains. In the mapping approach we 
lose information about the sequence of walls visited by the swimmer (though it can be deduced 
a posteriori), but the reduction in dimensionality is crucial to understanding the dynamics. For 
example, consider a microorganism inside a hexagonal billiard, with interior angle 9p = 120°, and 
departure angle 9c = 42°, so that the swimmer may arrive next at either the adjacent wall or 
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Figure 6: (a) For motion inside a hexagon with 9c £ (30°, 60°), a swimmer departing from one wall may collide 
either with the adjacent wall or the subsequent wall. A “beam” of all possible swimmers leaving from the same wall 
is shaded grey. Swimmers departing from x € (a, 1) collide with the adjacent wall, an unstable, stretching region; 
swimmers departing from x € (0, a) skip a wall and collide with a stable, focusing region, (b) The piecewise linear 
map for 9c = 42°. The slope indicates beam stretching or focusing. 


the subsequent wall, as illustrated in Fig. |^. The shaded region shows the space covered by the 
original beam of swimmers departing from the entire length of one wall. The map from [0,1] to 
itself (see Fig. |^-b) takes 0 to a point /?, and a point a to 0 (by convention due to the ambiguity 
at the vertex), and is given by 


Xn+l = f(Xn) = 


/3a ^(a - Xn) Xn < a, 

(1 - a)“^(l - Xn) Xn > a, 


(5) 


where (a,/3) ^ (0.54,0.37). The definition of /3 is just as in Eq. (lb), where it also corresponds to 


the image of the point x = 0. The values a and /3 for arbitrary (iV, 0c) are included in Appendix 


[a} Using the beam analogy, the region x G {a, 1) maps to the full length of the adjacent wall, 
corresponding to a decreased swimmer density upon arrival and a slope of the map with absolute 
value greater than one. Meanwhile, the region x G (0, a) maps to a region of length smaller than 
a, corresponding to an increase in the swimmer density and a slope of the map with absolute value 
less than one. 

The analogy illuminates individual swimmer dynamics: regions associated with a decreasing 
swimmer density, or stretching, are unstable for the swimmer trajectory, regions where the density 
of swimmers is unchanged are neutral, and regions of increasing density, or focusing, are stable. For 
example, the trajectory with departure angle 9c = 42° and initial point xq = 0.685 inside a hexagon 
is shown in Fig. The swimmer begins by reflecting seven times with walls adjacent to each 
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wall of departure. These reflections take place in the unstable region of the map, and the swimmer 
is eventually expelled to the other (stable) region where it rapidly settles to a consistently wall¬ 
skipping periodic triangular orbit. Note that this stable trajectory is precisely the one predicted 
for motion inside a regular triangle; the swimmer only visits three walls and cannot distinguish 
the hexagonal geometry from the triangular geometry. The stretch factor I3a~^ is less than one, 
resulting in stable periodic orbits. If the number of walls is prime, however, a stable wall-skipping 
orbit cannot trace out a regular polygon, and generally results in a more complex trajectory as 
observed in Fig. (solid line) for = 5 or Fig. for N = 7. 




Figure 7: (a) The trajectory of a swimmer with departure angle 6c = 42° and initial position xo = 0.685 in a 
hexagon {N = 6). After seven collisions with walls adjacent to each wall of departure the swimmer skips a wall and 
settles rapidly to a consistently wall-skipping periodic triangular orbit. The trajectory is colored red after the first 
30 reflections, (b) A stable, wall-skipping periodic orbit with a prime number of walls cannot trace out a regular 
polygon (here N = 7, 6c = 30°, and xq = x* = 0.093). 


Consider now a seemingly simpler example, microorganism billiards inside a square with 6c G 
(vr/d, 7r/2), as illustrated in Fig. [^. When the swimmer approaches the wall opposite the wall 
of departure it has changed its orientation relative to the mapping scheme, and the direction of 
increasing x must be reversed. The piecewise linear map for this case is given by 

{ 1 - (3a~^{a - Xn) Xn<a, 

( 6 ) 

(1 - a) (1 - Xn) Xn > a, 

which is now a continuous map due to the orientation reversal on the opposite wall (see Fig. [^). 
There remains one large unstable, stretching region corresponding to adjacent wall collisions, but 
now the opposite wall is a neutral region in which the density of swimmers in an incident beam 
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Figure 8: (a) Motion inside a square with 6c G (rr/d, 7r/2) is illustrated. Unlike in Fig. only neutral and unstable 
regions exist, and the resulting dynamics are chaotic, (b) An orientation reversal results in a continuous map from 
[0,1] to itself. 


is unchanged. The long-term dynamics of an individual swimmer are determined by the interplay 
between the unstable and neutral regions, and in particular the frequency of visits to the unstable 
stretching region and the degree of stretching there. Two such trajectories are shown with departure 
angles 9c = 52° and 9c = 72° in Fig.|^&:c. The dynamics are chaotic in both examples, but there is 
a clear distinction between the resulting trajectories. For 9c = 52°, after a few transient reflections 
(not shown) the trajectory is confined to a small region of the interior domain, while in the second 
case the dynamics do not at first glance appear to be contained in any such subregion. To study 
such long term dynamics in more detail we will proceed in the next section to explore the invariant 
measure and the Lyapunov exponent. 

For large, finite values of N , integer multiples of tt/N are special values of 9c for which the 
beam collides with only one other wall and all initial points xq give a neutrally stable periodic 
orbit. Figure 10 r shows such a finite estimate of such an example, with N = 200, 9c = 757r/200, 
and xq = 0.2. If the beam intersects two walls, let 9c = {a + b)~^ \a{kT: /N) + b{k + Ijvr/A^] for 
some a, 6 > 0, so that k = N9c/'x — b{a + b)~^. Then, for large N, 


{a,/3) 


a + b 




TTab 


(a -|- by tan{9c)N 


( 1 ,- 1 ), 


(7) 


so that 


I3a ^ ~ 1 — 


27ra 


[a -|- b) taii{9c)N 


< 1 . 


( 8 ) 
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Figure 9: (a) Simulation of a swimmer in a square with 6c = 52°. (b) Simulated prediction of the invariant measure 
for 6c = 52°. The final distribution shows that there are four domains that have nearly equal numbers of swimmers. 
Swimmers that start initially close together may end up in drastically different places with in these domains or 
possibly in entirely separate domains which makes this system chaotic, (c) Simulation of a swimmer in a square 
with 6c = 72°. (d) Simulated prediction of the invariant measure for 6c = 72°. The system appears nonuniform and 
depicts a seemingly ergodic system. 
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Figure 10: With N = 200 sides, and initial point xq = 0.2, periodic trajectories are found for the special angles: (a) 
8c = 757r/200 (period p = 8), and (b) 6c = 767r/200 (period p = 50). 


Hence, for large but finite N there remains a stable focusing region, but the period of the stable 
periodic orbit can be quite large. Figure [T^ shows a trajectory with N = 200, 6c = 767r/200, and 
xq = 0.2, where the period is found to be p = 50 (the smallest p for which 76p = 200m with m,p 
integers). 

Finally, in the limit as —>■ oo (a circular boundary) the entire boundary may be mapped 
onto itself via the internal angle, [0, 27r] —)• [0, 27r] by writing 0^+1 = 6n + 20c (assuming 6c < vr/2). 
Periodic trajectories are therefore found whenever 6c = {q/p)'^' for integer q and p, and then p 
is the period of the trajectory; otherwise, the map is chaotic but returns to within an arbitrarily 
small distance of any point on the circle infinitely often. 


2.3. Lyapunov exponent for the dynamics 

A useful means of characterizing the long-term behavior of the billiard system is to study the 
nature of the invariant measure for the dynamics, or the distribution of swimmers on each wall 
that is preserved with each iteration of the mapping. The invariant measure for the case of purely 


adjacent wall reflections (see 12.1) is simply the fixed point of the map, x* = f{x*) in Eq. (la). 
Writing the normalized invariant distribution as p{x)^ this case corresponds to a density with 
support at only one point, p{x) = 5x*{x), where 6x*{x) is a Dirac delta function centered at x*. 

The invariant measure for more irregular cases, such as for swimming in a square with 6c € 
(vr/d, 7r/2), is more complicated but speaks to the different possible structures observed in the 
long-time behavior of the system, for instance whether or not the trajectories are ergodic over the 
polygon’s interior (space-filling). The invariant measures for the cases 6c = 52° and 6c = 72°, 
depicted in Fig. [^,d, were computed by taking a uniform distribution of initial positions, xq ~ 
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U[0, 1], and measuring their locations after 4000 wall reflections (and finally, normalizing). In the 
first case, each individual swimmer is attracted to a small region of space, reflecting off of only 
four separate small domains. Depending on the initial position, the invariant set is one of four 
symmetric rotations of the region shown in Fig. [^. In the second case a single swimmer trajectory 
appears to cover the entire domain, as for a classical ergodic system. 

It is natural to ask about the rate at which two neighboring trajectories converge together or 
diverge from each other — the maximal Lyapunov exponent for the dynamical system. This is 
defined as 


A 


n—1 


lim - ^log|/'(xi)|. 


n^oo 77 , 


2=0 


(9) 


where Xj+i = f{xi). The distance along a wall between two swimmers after n reflections is 
approximately Axn ~ Axoexp(An), where Axq is their initial separation. For 9c = 52°, the 
Lyapunov exponent is approximately A = 0.12, while for 9c = 72° the value is larger, A = 0.47. From 
these two numbers, we can determine that the configuration in Fig. undergoes approximately a 
third as much stretching as that in Fig. [^. This tells us one of two things about the system with 
9c = 52°: either this system spends less time in the chaotic region or it has a chaotic region that 
does only a small amount of stretching per iteration. 




Figure 11: (a) In a square domain, the maximal Lyapunov exponent is exactly A = log(tan(0c)) for 9c < 7r/4, and 
then increases to a very noisy peak region before decaying again, always bounded above by 1. (b) The Lyapunov 
exponent is shown for trajectories internal to regular polygons with = 5,10 and 20 sides. Regions of overlap are 
departure angles for which the swimmer cannot distinguish between geometries. 


Figure [TT^ shows the computed Lyapunov exponent as a function of 9c for the square domain. 
For 9c < vr/4, the dynamics are confined to adjacent wall reflections, and from Eq. (la) we have that 


|/'(xj)| = tan(0c)- The stability of the dynamics in the adjacent wall case corresponds to negative 
values of A, and specifically A = log(/3) = log(tan(0c)) is exact for 9c < vr/4. For 9c G {'k j 
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however, the contribution to A depends on the frequency of reflections from the unstable region 
(where \ f{xi)\ > 1) and the neutral region (where \f'{xi)\ = 1). The Lyapunov exponent is very 
rough in the region where it is largest. 

Figure [T^ shows the computed Lyapunov exponents for three other polygonal domains, with 
= 5, 10, and 20. When 0^ < tt/N the swimmer only moves to the adjacent wall and A = log(/3) = 
log(sin(0c)/sin(27r/iV — 9c)) < 0, exactly. For one of the special angles 9c = kir/N with integer 
k and 9c < 7r/2, such that the departure wall maps perfectly onto one other wall, the periodic 
trajectory is neutrally stable and we have exactly A = 0. For other departure angles, the dynamics 
can be stable and periodic, or unstable and chaotic. Regions of overlap are departure angles for 
which the swimmer cannot distinguish between the geometries, as in the octagon/square example 
in Fig. I^and the hexagon/triangle example in Fig. For microorganism billiards inside a circle, 
the dynamics are always neutral, and the limiting behavior is A —)• 0 for all fixed 9c- 

Note that the computed value A is the Lyapunov exponent of the one-dimensional map, which 
differs from a “true” Lyapunov exponent since the transit time between reflections varies through¬ 
out the dynamics. Defining the nth transit time between reflections as tn, let the total time 
Tn = with ti < T. Here r is the maximum time between collisions, which depends on 

the diameter of the domain and the speed of the organism (r = diam/C/). The true Lyapunov 
exponent A (with units of inverse time) is 


. n—1 

A= lim Tf^^log\f'{xi)\, 
Tn ^ 


which differs from A (see Eq. Q) by the factor of rather than n in the denominator. We thus 
have 

Ti 

A = lim — A„, (10) 

n^oo 

where \n is Eq. without the limit. By Oseledec’s multiplicative ergodic theorem [^, both 
limits defining A and A must exist, which implies that the average time between encounters T = 
lim^^oo Tn/n < r must also exist. Hence, we have 


A = A/r, |A|>|A|/r. 


Since A and A have the same sign, chaos in the map is reflected as chaos in the full system, as is 
the absence of chaos. 


2-4- Robustness to random fluctuations 

The swimming trajectory of a microorganism is unlikely to be straight even over short distances. 
This can be due to internal biological mechanisms, such as the run-and-tumble flagellar dynamics of 
E. coli, thermal fluctuations, or other hydrodynamic effects. Alternatively, the angle of departure is 
likely to vary somewhat with each wall interaction, which was found to be true for Chlamydomonas 
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wall reflections by Kantsler et al. m- We would like to know whether the stable periodic orbits 
found for the adjacent wall case, with /3 < 1, survive the introduction of weak random fluctuations. 

As a hrst approximation to the effect of randomness we assume a Gaussian distribution around 
the deterministic point of arrival, writing for adjacent wall interactions in a regular polygon, 


Xn+l = /3(1 - Xn) + aZn+l, 


( 11 ) 


where j3 is given in Eq. (la), Zn ~ A^(0,1) are independent normal Gaussian random variables, 
and a is the standard deviation. The solution to © is 


Xn =(-/3)”'Xo 
n—1 

+ - aZ^_i) + aZn 

i=l 

=(-grxo + ;9 / 7 (12) 

^ i=0 

The last term is a sum of normal random variables with different variances, which leads to 


Xn 


{-/3Txo + I3 


1 - {-pr 

1 + /3 


+ cr 


/ 


1 - 

1-/32 


Zn, 


(13) 


where Zn ~ 3V(0,1). In the case /3 < 1, the transient dynamics decay exponentially fast and the 
trajectory is simply given by 


/? 


+ 


(7 Zn 


1 + /3 


n > 1. 


(14) 


The dynamics are still focused around the stable limit point in the deterministic case, but with 
a Gaussian spread that retains a memory of the recent past. When /3 ~ 1 the variance of the 
dynamics may be very large, but this is also the setting where the other walls should really be 
taken into consideration. Nevertheless, for reasonably small a and for adjacent wall reflections, the 
dynamics are robust, settling neatly into the deterministic periodic orbit with only small deviations. 

Another biologically relevant form of randomness, observed in the reflections of Chlamydomonas 
cells m, is a random departure angle. If the nth departure angle is given by 9c + (xZn, again with 
Zn ~ A^(0,1), then this results in a modified and random stretch factor j3 ^ j3n = sin(0c + 
aZn)l sin(0p + 6c + crZn)- While large tail events may lead to a non-adjacent wall reflection, which 
could result in a signihcant shift of the trajectory, we will neglect such events for the present 
calculation. The mapping is then given by Xn+i = /3n(l — Xn), resulting in a trajectory 


n 

Xn — QnXo ^ ^ Qj , 

1=1 


1-1 

Qj = (- 1 )'^ Pn-l-i, 


i=0 


(15) 

(16) 
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which may be used for computational purposes but is not analytically tractable. If, however, 
we assume a small variance and Taylor expand about small cr, then ~ /3 + where r = 

sin(0p)/sin^(0c + Op), and then inserting into the above and neglecting all terms of O(cr^) we find 




Qj = (-!)■’ + raZn-i-i) 


i=0 




V i=0 

= (- 1 )^ (^(3^ + P^~'^ray/jZj^ , 


n—l—i 


where Zj ~ A^(0,1). Then, after summing over normal variables as before, we have 


(_!)- (^/3-+ r/3"-X0 + 


/3(i-(-ir) 

1 + /3 


+ ra. 


'l - (1 + n(l - /32))/32n 


(l-/32)2 


where Z'^ ~ .^(0,1). For cr = 0 we recover the deterministic trajectory from Eq. Q, and for n 
large the trajectory is again drawn to the deterministic fixed point, with 


/3 1 

x„.^z -^ n>l. 


1 + /3 1 -/ 32 ' 


(17) 


Note that r(l — ^ = csc(0p + 20c), with Op = {N — 2 )tt/N the interior angle of the polygon. 


3. Passive sorting with boundary geometry 


With an eye towards the manipulation of microorganism populations for basic biological re¬ 
search and potential engineering applications |60j . numerous capture and sorting techniques have 
recently been designed using wedge-shaped boundaries (HU |62l E3] , chevron and heart-shaped chips 
[64| . smooth microchannels [HSl [Ml Wi\ . corrugated microchannels [M], non-convex boundary ge¬ 
ometries [69] , and regular patterned arrays [7011711 [Z2|. Using the hndings of the previous section 
we propose a new passive sorting technique for species with different departure angles. In par¬ 
ticular, we consider swimmers with departure angles 0c = 12° and 0c = 20°, which are the mean 
values for two different Chlamydomonas strains m- The sorting device is shown in Fig. a 


channel connects two square chambers of unit area which are relatively rotated by 90°, and for the 
lengths d and g shown in the hgure we choose d = 0.25 and g = 0.18. The centers of the squares 
are separated by a distance (3 -|- y/2)l2. Initially, 100 of each swimmer type are distributed with 
random initial position and orientation in the device, but at t = 10 (with swimmers moving at unit 
speed) the two strains are perfectly sorted. 
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t = 0 



Figure 12: Two hundred swimmers with unit speed, half with 6c = 12° (circles) and half with 9c = 20° (triangles), 
placed at random inside a simple sorting device. Their positions are shown at t = 0 (top), and then again perfectly 
sorted at t = 10 (bottom). 


The two swimmer types are separated by choosing the lengths d and g so that the fixed points 
of the map associated with each swimmer lie on a wall in one square and on a gap in the other. 
For a square domain, the fixed point is given by x*{9c) = tan(0c)/(l + tan(0c)) (see ^2.1| ). The 
two swimmers above are therefore successfully sorted when x*(12°) < d,g < x*{20°), or 0.175 < 
d,g < 0.267. The distance between the channels is not too important but without at least a small 
separation there can be more complicated effects near the opening in the leftmost square. If instead 
a rectangular channel (parallel walls) connects the square domains, it can be shown that the device 
sorts two swimmers so long as one has Be < 22.5° (exactly) and the other has 9c > 22.5°. 

Figure 13 shows the same sorting device but in the case that the departure angles vary with 
a Gaussian spread about the respective means with standard deviation a. The swimmer distri¬ 
butions are shown again at t = 10 in the cases with a = 2.5°, where the sorting is significant 
but imperfect, and with u = 5°, where the sorting is showing signs of deterioration. To measure 
the effectiveness of the sorting device as a function of the fluctuations in the departure angle, we 
consider swimmers with uniformly distributed initial positions and orientations inside the domain, 
and seek the probability that a swimmer is contained in the appropriate chamber after a long time. 
Specifically, we define an order parameter S = (Pi -I-P 2 ) — 1, where Pi is the probability that a cell 
with 0c = 12° is contained in the leftmost chamber as t —)• 00 and P 2 is the probability that a cell 
with 9c = 20° is contained in the rightmost chamber as t —>■ 00 . When S = 1 the system outside of 
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(7 = 2.5° 



Figure 13: The distribution of swimmers at t = 10 with Gaussian perturbations of the departure angle, with standard 
deviation (a) a = 2.5° and (b) (7 = 5°. 


a measure zero set is eventually sorted perfectly, while 5 = 0 is indicative of a disordered system. 
We estimate the order parameter by simulating 10, 000 randomly placed swimmers of each type 


and counting the number of properly sorted swimmers at t = 100. The result is shown in Fig. 14 


Intuitively we have perfect sorting for small Gaussian fluctuations of the departure angles, with a 
steadily diminishing sorting for large fluctuations. 


4. Microorganism billiards in a periodic array of square obstacles 

The approach used to study microorganism billiards inside a regular polygon can be fruitfully 
applied to the external problem of locomotion in a periodic array of polygonal boundaries. Recent 
related experiments by Volpe et al. m and Brown et al. [35] have shown ballistic, diffusive, and 
entrapped dynamics of synthetic swimming Janus particles in a lattice of obstacles. Battista et 
al. [73] studied the locomotion of malaria parasites, Plasmodium sporozoites, in an array of round 
pillars in a hexagonal lattice and connected stable migration to the organism’s mechanical flexibility. 
Here we consider a swimmer in an infinite lattice of square pegs, each of unit area with centers 
separated by a distance L = 1.65. The identification of each surface with every other surface 
once again results in a one-dimensional map that gives considerable insight into the swimming 
trajectories associated with a given departure angle 6c. As in the internal problem, Xn denotes 
the position of the swimmer’s nth reflection on any surface, where x is increasing in the direction 
of swimming. We consider only the case where the body does not slide along the surface before 
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Figure 14: The sorting order parameter as a function of the departure angle standard deviation. S = 1 indicates a 
perfectly sorted system, while S = 0 indicates a randomized swimmer distribution. 


departing, setting 5 = 0, and explore three cases that are modestly representative of the myriad 
possibilities in the external problem. 

Figure [TS^ shows the first representative example, a case where the departure angle 6c is such 
that a beam of swimmers departing from a given surface arrives on two surfaces of one other square 
obstacle. The one-dimensional map, which is continuous in this case, and a sample trajectory for 
6c = 25° are shown in Figs. [T^fcc. The map is composed of a focusing region (|/'(x)| < 1 for 
X £ [0, a]) and a neutral region {\f'{x)\ = 1 for x > a), where a = 0.26. The result is a trajectory 
that is at first ballistic while the body reflects from neutral surface to neutral surface, but eventually 
steps into the stable focusing region, corresponding here to a periodic trajectory. The resulting 
orbit is identical to the orbit in the internal problem in a square domain. Generally, if 6c < 45°, 
then the map can only contain neutral and stable regions, so that the microorganism billiard is 
eventually trapped into a periodic trajectory. In addition, reflections from parallel surfaces are 
neutral since they do not vary the length of the beam of swimmers arriving there, so any focusing 
or stretching dynamics must also be associated with directional changes (hitting walls that are 
perpendicular to the wall of departure). 

A second example is shown in Figs. & -c, where 6c = 48°. The swimmer can arrive at two 
surfaces, but now on two different obstacles. There remains a neutral region (|/^(x)| = 1 for 
X G [0,a]) with a = 0.41. However, since 6c > 45°, the region x G [a, 1] is now stretched 
(|/'(x)| > 1). The resulting trajectory therefore never settles to a stable periodic orbit. Figure [IB}; 
shows a sample trajectory for this case. Although there can be no stable periodic orbit, we do 
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Figure 15: The external microorganism billiard problem in an array of unit square obstacles, (a) Illustration of the 
case with obstacle distance L = 1.65 and departure angle 6c = 25°. The trajectory only reaches the two visible 
sides of one other obstacle when leaving from any given surface, (b) The one-dimensional map of the dynamics is 
continuous in this case, with a stable focusing region for x £ [0, a] and a neutral region for x £ [a, 1] for a = 0.26. (c) 
An illustrative trajectory. The swimmer is eventually trapped into a periodic orbit identical to that in the internal 
billiard problem in a square. 
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Figure 16: (a) Illustration of the case L = 1.65 and 6^ = 48°. (b) The map contains a neutral region for x £ [0, a] and 
an unstable stretching region for x G [a, 1], with a = 0.41. Reflections from parallel surfaces are neutral, so unstable 
stretching regions also correspond to directional changes, (c) An illustrative trajectory. The dynamics jump back 
and forth between the neutral and unstable regions of the map, resulting in a chaotic trajectory that nevertheless 
retains a coherent structure. 
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observe an order to the chaotic trajectory. It is possible, for instance, for the trajectory to undergo 
a periodic oscillation between the neutral and unstable regions of the map, leading to a weakly 
unstable periodic trajectory. 
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Figure 17: (a) Illustration of the case L = 1.65 and 6c = 62°. (b) The map contains two neutral regions corresponding 
to parallel wall reflections, and one unstable stretching region corresponding to directional changes, (c) An illustrative 
trajectory. The trajectory is chaotic; the swimmer is occasionally drawn into temporarily trapped dynamics, as in 
the Fig. [ISf, but due to the swimmer sampling the unstable stretching region. 


As a third and final example we set 9c = 62°, shown in Figs. dZh -c., with a = 0.65. In this case 
the swimmer can reach three surfaces; the two surfaces parallel to the initial wall are neutral, while 
the surface perpendicular to the initial wall is an unstable stretching region. As in the previous 
example, the trajectory therefore samples both the neutral and the stretching regions, leading 
to a chaotic trajectory. A quasi-ballistic drift along neutral (parallel) surfaces is punctuated by 
occasional directional changes (perpendicular surfaces). Counter-intuitively, the swimmer is thus 
occasionally drawn into what appear to be almost trapped dynamics, as in the first example, but 
this is due to the swimmer sampling the unstable stretching region. 
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5. Discussion 


In this paper we investigated the two-dimensional billiard-like motion of microorganisms which 
upon reflection from a surface depart with an angle independent of the angle of incidence. We 
first considered the swimming dynamics inside a regular N-sided polygon. For departure angles 
6c £ (0,7r/A^), resulting in only consecutive wall impacts, the swimmer was found to settle into a 
stable periodic trajectory (an inscribed regular polygon of the same degree). This stable orbit was 
shown to be robust to small random fluctuations of the arrival position or the departure angle, 
and independent of the sliding distance 6 so long as the body does not slide past a corner. What 
dynamics will ensue if the body does in fact slide past a corner remains unclear. Depending on the 
swimmer geometry and propulsive mechanism the swimming trajectory may be a simple extension 
of those already discussed, there may be trapping at the corner below a critical polygon interior 
angle, there may be continuous sliding along the entire domain [36], or possibly something more 
exotic will appear. Though not discussed here, trapping at corners with acute opening angles is 
possible even without sliding [38] . Even in classical billiards the effects of corners is an active topic 
of research [HES]. 

For non-consecutive impacts with 9c G vr/2), a one-dimensional piecewise linear map may 

be used to characterize the reflection regions as stable (focusing), neutral, or unstable (stretching), 
which is an easy way to predict whether the trajectory will settle to a stable orbit or undergo 
unstable chaotic dynamics. The Lyapunov exponent for these maps proved to be a useful measure 
of chaos since there is a direct relationship between the chaos in the map and the chaos in the full 
system. The special case in which the angle of departure is an integer multiple of the interior angle 
of the polygon {6c = kir/N) resulted in neutrally stable periodic orbits for any initial position. The 
stable fixed points derived for microorganism billiards inside a regular polygon were then used to 
design a model sorting device to passively separate two different types of swimmers. The effec¬ 
tiveness of the sorting device was investigated as a function of Gaussian variance in the departure 
angle. The present work may shed light on the nature of repeated wall reflections in conhned and 
patterned domains, and may suggest novel methods for directing the transport of microorganisms. 

The one-dimensional mapping approach was also applied to the external problem of a mi¬ 
croorganism swimming in an infinite array of square obstacles, for which we presented a few 
representative examples. Depending on the departure angle and the spacing between obstacles, 
the trajectory may become bound in a trapped orbit, or may undergo more complex and possibly 
chaotic dynamics. A more complete characterization of trajectories and diffusive behavior in the 
external problem remains an interesting future direction for study. 

To understand the behavior of larger microorganisms in confined domains will likely require 
more detailed modeling of the body mechanics and hydrodynamic interaction with the surface. 
For instance, experimental and numerical studies of swimming spermatozoa past an edge have 
shown that the scattering angle depends on a complex elastohydrodynamic interaction Ha ES]. 
Another example is the flexible malaria parasite, Plasmodium sporozoite, which is drawn with 
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preference to the obstacles that provide a better fit to its shape in a hexagonal array of round 
pillars m- Yet larger organisms such as the nematode C. elegans, swimming in a regular array 
of obstacles, also exhibit complex trajectories and gait changes depending on the obstacle spacing 
While such organisms may be hydrodynamically attracted to surfaces and drawn inexorably 
towards one of infinite expanse IZSIE], they may slide off of the obstacles with sufficiently large gap 
spacing in the external problem, and therefore might still be well understood using the mapping 
approach presented herein. Looking forward, applications of microorganism billiard dynamics 
may therefore include improved spermatozoan selection to increase the success rates in in-vitro 
fertilization techniques, and the filtration of certain parasites for study or elimination by their 
mechanical properties. 

The authors thank David Anderson, Dwight Barkley, Renske Gelderloos, Michael Graham, 
Thomas Kurtz, and Douglas Weibel for helpful discussions. This research was supported by NSF 
grants DMS-1109315 and DMS 1147523. 

Appendix A. Exact expressions for the microorganism billiard inside a regnlar poly¬ 
gon 

For a regular polygon of degree N, treating the plane in complex variables and the side of 
departure emanating from the origin in the direction cq = 1, the sides are parallel to the unit 
vectors 


_ ^27vim/N 


(A.l) 


Assume that 9c G {kn/N, (k + Ijvr/A) for some integer k. Then, writing the swimming direction 
as dc = exp{i9c), we have 

k k 

^ ^ '^m T (dVk+1 — ^idcj ^ ^ Vm — O; T ^2dci 

m=0 m=0 


for as yet unknown transit distances Ai,A 2 . Multiplying each side in both equations by dc and 
taking the imaginary parts, we find 


which, using 


idc ^ ^2nim/N _|_ p^2'Ki{k+l)/_ q 

= —asin(0c), 


^m=0 


^ ^ ^'KirajN 

m=0 
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eventually simplifies to give 


where 


a = K 


P = n— 


sin {9c — kn/N) 

sm{9c) 

sin {9c — k'a/N) 
sin {2{k + 1)t:/N — 9c) 


sin {{k + l)7r/iV) 
sin {tt/N) 


Note that the important quantity fia ^ simplifies to 
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